Python数值插值和曲线拟合

亚图跨际

已于 2022-04-27 14:54:31 修改

阅读量2.6k

点赞数

分类专栏： Python 数值处理文章标签： python 数值插值和拟合

于 2021-03-28 06:00:01 首次发布

本文链接：https://blog.csdn.net/jiyotin/article/details/115275295

版权

Python 同时被 2 个专栏收录

400 篇文章 25 订阅

订阅专栏

数值处理

12 篇文章 1 订阅

订阅专栏

通常参与技术计算。数据来源可以是实验观察或数值计算。插值和曲线拟合之间是有区别的。在插值中，我们通过数据点构造一条曲线。这样做时，我们隐含了一个假设，即数据点是准确且不同的。相反，曲线拟合应用于通常由测量误差引起的包含散射（噪声）的数据。在这里，我们要找到一条平滑的曲线，可以从某种意义上近似该数据。因此，曲线不一定会击中数据点。插值和曲线拟合之间的差异如图1所示。

图1

多项式插值

拉格朗日法

插值的最简单形式是多项式。始终可以构造一个唯一的次数为 $n$ 的多项式，该多项式通过 $n + 1$ 个不同的数据点。获得该多项式的一种方法是拉格朗日公式，

$P_{n}(x)=\sum_{i=0}^{n} y_{i} \ell_{i}(x),\qquad(1a)$

其中下标 $n$ 表示多项式的次数，而

$\begin{aligned}\ell_{i}(x) &=\frac{x-x_{0}}{x_{i}-x_{0}} \cdot \frac{x-x_{1}}{x_{i}-x_{1}} \cdots \cdots \frac{x-x_{i-1}}{x_{i}-x_{i-1}} \cdot \frac{x-x_{i+1}}{x_{i}-x_{i+1}} \cdots \frac{x-x_{n}}{x_{i}-x_{n}} \\&=\prod_{j=0 \atop j \neq i}^{n} \frac{x-x_{i}}{x_{i}-x_{j}}, \quad i=0,1, \ldots, n\end{aligned},\qquad(1b)$

被称为基数函数。

例如，如果 $n + 1$ ，则插值是直线 $P_{1}(x)=y_{0} \ell_{0}(x)+y_{1} \ell_{1}(x)$ ，其中

$\ell_{0}(x)=\frac{x-x_{1}}{x_{0}-x_{1}} \quad \ell_{1}(x)=\frac{x-x_{0}}{x_{1}-x_{0}}$

当 $n = 2$ 时，插值是抛物线形： $P_{2}(x)=y_{0} \ell_{0}(x)+y_{1} \ell_{1}(x)+y_{2} \ell_{2}(x)$ ，其中

$\begin{array}{l}\ell_{0}(x)=\frac{\left(x-x_{1}\right)\left(x-x_{2}\right)}{\left(x_{0}-x_{1}\right)\left(x_{0}-x_{2}\right)} \\\ell_{1}(x)=\frac{\left(x-x_{0}\right)\left(x-x_{2}\right)}{\left(x_{1}-x_{0}\right)\left(x_{1}-x_{2}\right)} \\\ell_{2}(x)=\frac{\left(x-x_{0}\right)\left(x-x_{1}\right)}{\left(x_{2}-x_{0}\right)\left(x_{2}-x_{1}\right)}\end{array}$

基数函数是次数为n的多项式，并具有以下性质

$\ell_{i}\left(x_{j}\right)=\left\{\begin{array}{l}0 \text { if } i \neq j \\1 \text { if } i=j\end{array}\right\}=\delta_{i j},\qquad(2)$

其中 $\delta_{i j}$ 是Kronecker增量。图2中针对 $x_{0}=0, x_{1}=2$ 和 $x_{2}=3$ 的三点插值（ $n = 2$ ）演示了此属性。

图2

为了证明插值多项式通过数据点，我们将 $x=x_{j}$ 代入等式（1a），然后使用等式（2）。结果是

$P_{n}\left(x_{j}\right)=\sum_{i=0}^{n} y_{i} \ell_{i}\left(x_{j}\right)=\sum_{i=0}^{n} y_{i} \delta_{i j}=y_{j}$

可以证明多项式插值的误差为

$f(x)-P_{n}(x)=\frac{\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n}\right)}{(n+1) !} f^{(n+1)}(\xi),\qquad(3)$

$\xi$ 位于区间 $\left(x_{0}, x_{n}\right)$ 的某处；其值是未知的。请注意，数据点离 $x$ 越远，则它对 $x$ 处的误差的贡献就越大。

牛顿法

尽管拉格朗日法从概念上讲很简单，但它并不适合有效的算法。使用牛顿法可以获得更好的计算过程，其中插值多项式以以下形式编写：

$P_{n}(x)=a_{0}+\left(x-x_{0}\right) a_{1}+\left(x-x_{0}\right)\left(x-x_{1}\right) a_{2}+\cdots+\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{n-1}\right) a_{n}$

该多项式有助于进行有效的评估程序。例如，考虑四个数据点（ $n = 3$ ）。这里的插值多项式是

$\begin{aligned}P_{3}(x) &=a_{0}+\left(x-x_{0}\right) a_{1}+\left(x-x_{0}\right)\left(x-x_{1}\right) a_{2}+\left(x-x_{0}\right)\left(x-x_{1}\right)\left(x-x_{2}\right) a_{3} \\&=a_{0}+\left(x-x_{0}\right)\left\{a_{1}+\left(x-x_{1}\right)\left[a_{2}+\left(x-x_{2}\right) a_{3}\right]\right\}\end{aligned}$

可以通过以下递归关系向后求值：

$\begin{array}{l}P_{0}(x)=a_{3} \\P_{1}(x)=a_{2}+\left(x-x_{2}\right) P_{0}(x) \\P_{2}(x)=a_{1}+\left(x-x_{1}\right) P_{1}(x) \\P_{3}(x)=a_{0}+\left(x-x_{0}\right) P_{2}(x)\end{array}$

对于任意 $n$ 我们有

$P_{0}(x)=a_{n} \quad P_{k}(x)=a_{n-k}+\left(x-x_{n-k}\right) P_{k-1}(x), \quad k=1,2, \ldots, n,\qquad(4)$

用xData表示数据点的 $x$ 坐标数组，用 $n$ 表示多项式的次数，我们有以下算法来计算 $P_{n}(x)$ ：

p = a[n]
	for k in range(1,n+1):
		p = a[n-k] + (x - xData[n-k])*p

通过强制多项式通过每个数据点来确定 $P_{n}$ 的系数： $y_{i}=P_{n}\left(x_{i}\right), \quad i=0,1, \ldots, n$ 。这产生了这些联立方程：

$\begin{aligned}&\begin{array}{l}y_{0}=a_{0} \\y_{1}=a_{0}+\left(x_{1}-x_{0}\right) a_{1} \\y_{2}=a_{0}+\left(x_{2}-x_{0}\right) a_{1}+\left(x_{2}-x_{0}\right)\left(x_{2}-x_{1}\right) a_{2}\end{array}\\& \vdots\\& y_{n}=a_{0}+\left(x_{n}-x_{0}\right) a_{1}+\cdots+\left(x_{n}-x_{0}\right)\left(x_{n}-x_{1}\right) \cdots\left(x_{n}-x_{n-1}\right) a_{n}\end{aligned},\qquad(a)$

引入除数差

$\begin{aligned}&\begin{aligned}\nabla y_{i} &=\frac{y_{i}-y_{0}}{x_{i}-x_{0}}, \quad i=1,2, \ldots, n \\\nabla^{2} y_{i} &=\frac{\nabla y_{i}-\nabla y_{1}}{x_{i}-x_{1}}, \quad i=2,3, \ldots, n \\\nabla^{3} y_{i} &=\frac{\nabla^{2} y_{i}-\nabla^{2} y_{2}}{x_{i}-x_{2}}, \quad i=3,4, \ldots n\end{aligned}\\&\qquad\qquad\vdots\\&\nabla^{n} y_{n}=\frac{\nabla^{n-1} y_{n}-\nabla^{n-1} y_{n-1}}{x_{n}-x_{n-1}}\end{aligned},\qquad(5)$

等式（a）的解为

$a_{0}=y_{0} \quad a_{1}=\nabla y_{1} \quad a_{2}=\nabla^{2} y_{2} \quad \cdots \quad a_{n}=\nabla^{n} y_{n},\qquad(6)$

如果系数是手工计算的，则使用下表1中的格式（ $n = 4$ 所示）很方便。

表1

表中的对角项（ $y_{0}, \nabla y_{1}, \nabla^{2} y_{2}, \nabla^{3} y_{3}$ 和 $\nabla^{4} y_{4}$ ）是多项式的系数。如果以不同的顺序列出数据点，则表中的条目将发生变化，但所得多项式将是相同的-请记住， $n$ 次插值 $n + 1$ 个不同数据点的多项式是唯一的。

可以使用以下算法在一维数组 $a$ 中进行机器计算（我们使用符号 $m = n + 1 =$ 数据点数）：

a = yData.copy()
for k in range(1,m):
		for i in range(k,m):
				a[i] = (a[i] - a[k-1])/(xData[i] - xData[k-1])

最初，a包含数据的y坐标，因此它与上表1中的第二列相同。每次通过外部循环都会在下一列中生成条目，这些条目将覆盖a的相应元素。因此，最终包含表1的对角项（即多项式的系数）。

牛顿插值算法代码

newtonPoly模块包含牛顿法进行插值所需的两个函数。给定数据点数组xData和yData，函数coeffts返回系数数组a。找到系数后，可以使用函数evalPoly在x的任何值上评估内插 $P_{n}(x)$ 。

内维尔法

牛顿的插值方法包括两个步骤：系数的计算，然后是多项式的求值。如果使用相同的多项式在 $x$ 的不同值处重复执行插值，则效果很好。如果仅要插值一点，则使用内维尔算法一步计算插值的方法是更方便的选择。

令 $P_{k}\left[x_{i}, x_{i+1}, \ldots, x_{i+k}\right]$ 表示 $k$ 次多项式通过 $k + 1$ 个数据点 $\left(x_{i}, y_{i}\right),\left(x_{i+1}, y_{i+1}\right), \ldots,\left(x_{i+k}, y_{i+k}\right)$ 。对于单个数据点，我们有

$P_{0}\left[x_{i}\right]=y_{i},\qquad(7)$

基于两个数据点的插值是

$P_{1}\left[x_{i}, x_{i+1}\right]=\frac{\left(x-x_{i+1}\right) P_{0}\left[x_{i}\right]+\left(x_{i}-x\right) P_{0}\left[x_{i+1}\right]}{x_{i}-x_{i+1}}$

很容易验证 $P_{1}\left[x_{i}, x_{i+1}\right]$ 通过两个数据点; 也就是说，当 $x=x_{i}$ 时， $P_{1}\left[x_{i}, x_{i+1}\right]=y_{i}$ ，当 $x=x_{i+1}$ 时， $P_{1}\left[x_{i}, x_{i+1}\right]=y_{i+1}$ 。

三点插值是

$P_{2}\left[x_{i}, x_{i+1}, x_{i+2}\right]=\frac{\left(x-x_{i+2}\right) P_{1}\left[x_{i}, x_{i+1}\right]+\left(x_{i}-x\right) P_{1}\left[x_{i+1}, x_{i+2}\right]}{x_{i}-x_{i+2}}$

为了证明该插值确实与数据点相交，我们首先将 $x=x_{i}$ 代入，得到

$P_{2}\left[x_{i}, x_{i+1}, x_{i+2}\right]=P_{1}\left[x_{i}, x_{i+1}\right]=y_{i}$

同样， $x=x_{i+2}$ 得出

$P_{2}\left[x_{i}, x_{i+1}, x_{i+2}\right]=P_{2}\left[x_{i+1}, x_{i+2}\right]=y_{i+2}$

最后，当 $x=x_{i+1}$ 时，我们有

$P_{1}\left[x_{i}, x_{i+1}\right]=P_{1}\left[x_{i+1}, x_{i+2}\right]=y_{i+1}$

于是

$P_{2}\left[x_{i}, x_{i+1}, x_{i+2}\right]=\frac{\left(x_{i+1}-x_{i+2}\right) y_{i+1}+\left(x_{i}-x_{i+1}\right) y_{i+1}}{x_{i}-x_{i+2}}=y_{i+1}$

建立了模式之后，我们现在可以推导出通用的递归公式：

$\begin{aligned}& P_{k}\left[x_{i}, x_{i+1}, \ldots, x_{i+k}\right]=\\& \frac{\left(x-x_{i+k}\right) P_{k-1}\left[x_{i}, x_{i+1}, \ldots, x_{i+k-1}\right]+\left(x_{i}-x\right) P_{k-1}\left[x_{i+1}, x_{i+2}, \ldots, x_{i+k}\right]}{x_{i}-x_{i+k}}\end{aligned},\qquad(8)$

给定 $x$ 的值，可以按照以下表格格式（针对四个数据点显示）进行计算：

表2

用 $m$ 表示数据点的数量，计算表元素的算法为

该算法适用于一维数组y，该数组最初包含数据的y值（表2中的第二列）。每次通过外部循环，都会计算下一列中的y元素，这些元素将覆盖先前的条目。在过程的最后，y包含表格的对角线项。穿过所有数据点的插值的值（在x处评估）是y的第一个元素

内维尔插值算法代码

多项式插值的局限性

多项式插值应使用最少可行的数据点数进行。如果数据点之间的距离很近，则使用最接近的两个点进行线性插值通常就足够了。在大多数情况下，三到六个最近的邻居点会产生良好的结果。相交超过六个点的插值器必须令人怀疑。原因是远离兴趣点的数据点不会影响插值的准确性。实际上，它们可能是有害的。

使用过多点的危险如图3所示。圆圈表示11个等距的数据点。实线是内插值，它是所有点相交的次数为10的多项式。如图所示，如此高的多项式有在数据点之间过度振荡的趋势。通过使用跨越四个最近邻点的三次插值，可以获得更加平滑的结果。

图3

多项式外推（在数据点范围之外进行内插）是危险的。作为示例，请考虑图4。六个数据点显示为圆圈。第五级插值多项式由实线表示。插值在数据点的范围内看起来很好，但是当x> 12时，明显偏离明显的趋势。例如，在这种情况下，将x推算为14时，y是不合理的。

图4

如果无法避免外推，则可以使用以下措施：

绘制数据并直观地验证外推值是否有意义。
使用基于最近邻居数据点的低阶多项式。例如，对于图4中的数据，线性或二次插值将产生 $y (14)$ 的合理估计。
处理对数x与对数y的关系图，该图通常比x-y曲线平滑得多，因此推断起来更安全。通常，该图几乎是一条直线。这在图5中进行了说明，该图表示图4中数据的对数图。

图5

示例1

给定数据点

使用拉格朗日法确定 $x = 1$ 处的 $y$

$\begin{array}{l}\ell_{0}=\frac{\left(x-x_{1}\right)\left(x-x_{2}\right)}{\left(x_{0}-x_{1}\right)\left(x_{0}-x_{2}\right)}=\frac{(1-2)(1-3)}{(0-2)(0-3)}=\frac{1}{3} \\\\\ell_{1}=\frac{\left(x-x_{0}\right)\left(x-x_{2}\right)}{\left(x_{1}-x_{0}\right)\left(x_{1}-x_{2}\right)}=\frac{(1-0)(1-3)}{(2-0)(2-3)}=1 \\\\\ell_{2}=\frac{\left(x-x_{0}\right)\left(x-x_{1}\right)}{\left(x_{2}-x_{0}\right)\left(x_{2}-x_{1}\right)}=\frac{(1-0)(1-2)}{(3-0)(3-2)}=-\frac{1}{3} \\\\y=y_{0} \ell_{0}+y_{1} \ell_{1}+y_{2} \ell_{2}=\frac{7}{3}+11-\frac{28}{3}=4\end{array}$

示例2

数据点

在多项式上。通过构造一个类似于表1的除差表来确定该多项式的次数。

这是用于得出表中数字的一些示例计算：

$\begin{aligned}\nabla y_{2} &=\frac{y_{2}-y_{0}}{x_{2}-x_{0}}=\frac{59-(-1)}{4-(-2)}=10 \\\\\nabla^{2} y_{2} &=\frac{\nabla y_{2}-\nabla y_{1}}{x_{2}-x_{1}}=\frac{10-1}{4-1}=3 \\\\\nabla^{3} y_{5} &=\frac{\nabla^{2} y_{5}-\nabla^{2} y_{2}}{x_{5}-x_{2}}=\frac{-5-3}{-4-4}=1\end{aligned}$

从表中我们可以看到，牛顿多项式的最后一个非零系数（最后一个非零对角线项）为 $\nabla^{3} y_{3}$ ，这是三次项的系数。因此，多项式是三次的。

示例3

给定数据点

用内维尔法确定 $y (x) = 0$ 的根。

这是逆插值的示例，其中 $x$ 和 $y$ 的角色互换。我们不是在给定的 $x$ 上计算 $y$ ，而是找到与给定 $y$ 对应的 $x$ （在这种情况下， $y = 0$ ）。使用表2的格式（当然 $x$ 和 $y$ 互换了），我们得到

以下是表中使用的示例计算：

表格中的所有 $P$ 都是根的估计值，根是由涉及不同数据点的不同插值顺序得出的。例如， $P_{1}\left[y_{0}, y_{1}\right]$ 是根据前两个点通过线性插值获得的根，而 $P_{2}\left[y_{1}, y_{2}, y_{3}\right]$ 是使用后三个点进行二次插值的结果。从所有四个数据点的三次插值获得的根是 $x=P_{3}\left[y_{0}, y_{1}, y_{2}, y_{3}\right]=3.8317$ 。